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Managementuittreksel 

TNO  1  Kennis  voor  zaken 

SAS  interferometrie  voor  detectie  en 
classificatie 


Interferometrische  SAS  processing  is  ontwikkeld,  getest  en  toegepast  op 
data  uit  simulaties  en  uit  een  experiment  met  een  sonar  op  een  rail 
De  resultaten  laten  realistische  hoogtes  van  de  testobjecten  zien. 
Als  toepassing  wordt  detectie  van  begraven  objecten  beoogd,  maar  de 
methode  is  ook  geschikt  als  hulpmiddel  bij  classificatie  van  bodemobjecten 
en  voor  bathymetrische  survey. 


Probleemstelling 

De  akoestische  detectie  van  in  de  zeebodem 
begraven  objecten  vormt  een  uitdaging, 

/.elfs  voor  modeme  sonarsystemen  lnter- 
ferometnsche  synthetische  apertuur  sonar 
(inSAS)  wordt  gezien  als  een  veelbelovende 
techniek  om  deze  uitdaging  het  hoofd  te 
bieden.  Met  behulp  van  inSAS  kan  onder- 
scheid  gemaakt  worden  tussen  reflecties 
van  verschillende  hoogtes  en  vvellicht  kan 
de  bodemecho  onderdrukt  worden  ten 
gunste  van  echo’s  van  begraven  objecten 
Zodoende  kan  naar  vervvachting  de 
detectiekans  verhoogd  worden 


In  het  project  Actieve  SAS  technology  is 
daarom  gekozen  voor  het  ontvvikkelen  en 
implementeren  van  inSAS  Het  vverk  bouvvt 
voort  op  de  in  voorgaande  jaren  ontwikkel- 
de  softwaresuite  voor  SAS  processing,  be- 
wegingscorrectie  en  simulatie  Het  rapport 
geeft  de  status  van  de  ontvvikkehngen  vveer 
als  tussenrapportage 
Het  project  maakt  deel  uit  van  het  doel- 
financieringsprogramma  V512,  Sonar  en 
onderwater  propagatie  en  de  opdrachtgever 
is  de  Defensie  Mateneel  Organisatie 


Beschrijving  van  de 
werkzaamheden 

Om  interferometrie  toe  te  passen  /ijn  data 
benodigd  van  twee  parallelle  sonararrays  op 
verschillende  hoogte  die  de/el fde  sc£ne 
afbeelden.  Uit  faseverschillen  tussen  beide 
afbeeldingen  kunnen  hoogteverschillen  in 
de  scene  worden  afgeleid 
Inter ferometnsche  processing  verloopt  via 
de  tussenstappen  migratie  de  beeldvorm- 
ende  techniek  die  voor  beide  arrays  een 
fase-consi stent  beeld  moet  ople\  eren;  co- 
registratie,  waarbij  de  afbeeldingen  /o 
nauwkeung  mogehjk  over  elkaar  heen 
worden  gepast;  interferogram  berekenen. 
waarbij  de  afbeeldingen  ‘van  elkaar 
afgetrokken’  worden,  faseherstel  {phase 
unwrapping),  waarbij  fasesprongen  van 
meer  dan  twee  pi  ongedaan  gemaakt 
w  orden  en  ten  slotte  het  berekenen  van  een 
hoogtekaart’,  waarbij  de  faseverschillen  op 
basis  van  geometrische  gegevens  worden 
omgerekend  naar  hoogteverschillen  Voor 
elk  van  deze  fasen  is  software  ontwikkeld 
en  beproefd  Hiervoor  z.ijn  zowel  synthe- 
tische  data  uit  simulaties  gebruikt  als  expe- 
rimentele  data  van  metingen  op  een  rail 

Resultaten  en  conclusies 

Het  resultaat  van  het  beschreven  vverk  is 
een  inSAS  processing  suite  vvaarmee  rela¬ 
tive  hoogtes  in  een  sonarbeeld  berekend 
kunnen  worden  De  berekenmgen  met 
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hoogfrequente  data  geven  objecthoogtes  die 
realistisch  zijn  De/e  resultaten  geven 
vertrouvven  in  de  processing.  Op  diverse 
punten  is  de  processing  nog  voor  verbete- 
ring  vatbaar,  bijvoorbeeld  robuustheid  in 
aamve/igheid  \an  multipad-propagatie  en 
faseconsistentie  van  de  SAS-beelden 
Ook  de  toepassing  op  laagfrequente  data  /al 
nog  nader  onderzocht  worden 

Toepasbaarheid 

In  het  algemeen  is  de  interferometrische 
processing  ontvvikkeld  om  hoogteverschil- 
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Summary 


This  report  presents  an  overview  of  the  theory  and  implementation  of  interferometric 
S  AS  processing  at  TNO.  The  objective  of  the  technique  is  mainly  to  provide  a  possible 
solution  for  the  detection  of  buried  targets  in  mine-hunting  or  harbour  protection 
scenarios,  by  distinguishing  between  echoes  from  objects  at  different  heights,  such  as 
the  sea  bottom  surface  and  buried  objects  within  that  bottom.  Other  applications  are 
survey,  and  measurement  of  bathymetry'  in  general.  In  mine-hunting,  height  information 
from  interferometry  may  improve  the  detection  and  classification  of  naval  mines. 

The  processing  steps  in  interferometry,  required  to  derive  the  relative  heights,  are  the 
following: 

•  Generate  a  phase-consistent  (SAS-)image  for  each  receive  array.  Not  all  S AS 
processing  techniques  are  phase-preserving,  so  an  adequate  technique,  including 
motion  compensation  /  autofocusing,  should  be  selected. 

•  Co-register  (align)  the  images,  based  on  information  about  the  geometry'  and  on 
correlation.  Geometry -based  correction  is  usually  applied  to  the  whole  image, 
whereas  correlation  based  co-registration  is  often  applied  to  patches  of  the  image. 
The  dimensions  of  these  patches  have  to  be  chosen  carefully. 

•  Calculate  the  interferogram  by  cross-correlation  of  the  two  images, 

•  Undo  possible  phase  jumps  due  to  two  pi  ambiguities  to  obtain  absolute  phases, 
also  known  as  phase  unwrapping.  In  many  applications,  depending  on  the  height 
variations  and  geometry,  this  is  not  a  crucial  step.  Many  two-dimensional  phase 
unwrapping  methods  exist,  with  differences  in  sensitivity  to  noise,  robustness  and 
calculation  time. 

•  Calculate  the  height  map  using  geometrical  information  to  convert  absolute  phase 
differences  to  height  differences. 

For  each  of  these  steps,  processing  softw  are  has  been  developed  and  tested.  Testing  has 
been  carried  out  for  both  synthetic  data  from  simulations  and  measured  data  from  a  rail 
experiment. 

Promising  results  have  been  achieved  with  the  developed  processing.  Both  for 
simulated  and  for  measured  data,  correct  object  heights  are  found,  supporting  the 
validity  of  the  processing.  Remaining  improvements  lie  mainly  in  motion  compensation 
and  multipath  mitigation,  which  both  are  aimed  at  improving  the  phase  consistency. 

The  interferometric  processing  can  be  used  for  bathymetric  mapping  and  as  an  aid  in 
target  classification.  The  ultimate  objective  is  to  improve  the  detection  and 
classification  of  buried  mines,  although  this  concept  will  have  to  be  demonstrated. 
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1  Introduction 


Interferometry  is  a  technique  that  is  widely  applied  in  scientific  measurements  to 
measure  small  differences  in  height  or  length.  Optical  applications  are  well-established. 
In  sonar,  and  also  in  radar  (SAR,  synthetic  aperture  radar),  interferometry  is  mainly 
used  to  determine  height  differences  in  a  scene  under  observation  [1],  [2].  SAR 
applications  often  use  multiple  passes  of  the  radar  platform  to  determine  small  changes 
with  time  in  the  surveyed  area.  These  changes  may,  for  instance,  be  due  to  various 
geophysical  processes,  such  as  volcanic  activity  or  water  pressure  on  dams. 

This  is  known  as  repeat  pass  interferometry  and  can  be  achieved  with  only  one  receive 
antenna.  Other  applications  of  radar  interferometry  exist,  but  these  are  outside  the  scope 
of  this  report.  Owing  to  the  similarities  between  radar  and  sonar,  techniques  from  radar 
can  often  be  applied  to  sonar  as  well. 

Sonar  applications  of  interferometry  are  currently  mainly  found  in  bathymetric  surveys, 
where  variations  in  sea  bottom  topography  (bathymetry)  are  mapped  using  interfero¬ 
metric  side  scan  sonar,  based  on  single-beam  processing.  This  type  of  interferometry 
requires  at  least  two  arrays,  which  are  vertically  separated  by  a  fixed  distance,  referred 
to  as  the  baseline.  SAS  interferometry  is  quite  similar  as  it  is  also  applied  to  data  from 
two  receive  arrays.  The  difference  is  that  this  processing  is  applied  to  images  rather  than 
single  beams.  Height  differences  in  the  scene  under  observation  result  in  small 
differences  between  the  SAS-images  from  each  array.  Using  information  about  the 
geometry,  these  differences  can  be  translated  to  height  relative  to  a  plane  of 
Reference  [3],  [4],  [5],  [6], 

Current  international  developments  in  interferometric  SAS  are  mainly  aimed  at 
improving  mine  warfare  capabilities,  where  height  provides  additional  information 
about  detected  targets  that  helps  classification.  In  particular  on  autonomous  platforms 
(AUV’s:  autonomous  underwater  vehicles)  [7],  [8],  automatic  classification  is  essential 
in  achieving  the  desired  level  of  autonomy.  Hence,  all  available  extra  information  on 
detected  targets  is  useful.  However,  these  are  applications  where  high-frequency,  high- 
resolution  sonar  is  used  for  detection  of  proud  mines.  An  important  shortfall  in  mine 
warfare  capabilities  is  still  the  detection  of  buried  mines,  which  remain  hidden  in  high 
frequency  sonar  images.  Detection  in  sediment  (sand,  mud)  requires  lower  frequencies 
(roughly  between  5  and  30  kHz,  depending  on  sediment  type),  but  at  these  frequencies 
resolution  is  much  worse  as  a  consequence  of  smaller  bandwidth.  Low-frequency  SAS 
can  then  be  applied  to  improve  along-track  resolution  and  signal-to-noise-ratio  (SNR). 
In  a  single  SAS-image,  no  distinction  can  be  made  between  objects  at  the  bottom 
surface  or  below  the  surface.  Interferometry  then  may  provide  the  decisive  information 
whether  or  not  the  detected  object  is  buried  (i.e.  it  has  negative  height). 

Using  the  sea  floor  as  the  reference  height,  a  further  possibility  is  to  try  and  suppress  all 
scattering  from  the  bottom  interface,  for  example  by  subtracting  one  image  from  the 
other.  If  this  can  be  properly  achieved  (depending  on  bottom  bathmetry  and  resolution, 
among  others),  only  signals  from  below  or  above  the  bottom  will  remain,  thereby 
enhancing  the  echo  from  buried  targets.  This  TNO-patented  approach  is  known  as  the 
SUBDIRECT  algorithm,  which  can  be  investigated  once  the  interferometric  processing 
is  fully  functional  To  test  this  approach  experimentally,  ‘conventional1  interferometry 
has  to  be  developed  first. 
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In  the  TNO-project  ‘Active  SAS  Technology’,  the  development  of  interferometric  SAS 
(inSAS)  has  been  started.  The  goal  is  to  develop  inSAS  processing  for  the  detection  of 
buried  targets,  as  explained  above.  This  report  presents  the  current  status  of  the 
development.  Chapter  2  contains  an  overview  of  the  theoretical  background  and 
required  processing  steps  in  SAS  interferometry.  The  implementation  of  the  theory  in 
softw  are  has  been  tested  on  two  types  of  data,  simulated  and  measured. 

Chapter  3  presents  results  obtained  with  simulated  data;  Chapter  4  the  results  with 
measured  data.  The  measured  data  have  been  obtained  from  a  rail  experiment  in  1999 
by  DERA  (currently  QinetiQ)  and  GESMA,  from  the  UK  and  France,  respectively. 

Both  the  experimental  set-up  and  the  results  are  presented  in  Chapter  4. 

The  report  concludes  with  Chapter  5,  w  hich  gives  a  review  of  the  results  and  an  outlook 
to  further  developments. 
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2  Interferometry 


2.1  The  concept  of  interferometry 

The  objectives  of  SAS  interferometry  are  to  obtain  information  on  the  height  of  the 
seafloor  and  to  improve  the  classification  of  objects  by  retrieving  information  on  their 
height  and  -if  possible-  shape  (the  latter  depending  on  target  dimensions  and 
orientation,  scattering  strength,  etc.).  Interferometry  requires  that  data  are  acquired  by  at 
least  two  arrays  (Figure  2.1).  These  are  separated  by  the  baseline  vector  B.  As  shown  in 
Figure  2.1,  this  recording  geometry  enables  us  to  detect  subtle  time  differences  between 
arrivals  on  array  1  and  array  2  that  are  caused  by  elevation  differences. 


/TT 

Figure  2.1  Configuration  for  SAS  inlerferometry. 


In  the  following,  we  assume  that  only  one  transmit  element  is  used,  i.e  phase 
differences  only  occur  at  the  receiver  side.  We  refer  to  the  image  constructed  from  the 
data  acquired  by  array  1  as  the  master  image;  the  image  corresponding  to  array  2  is 
referred  to  as  the  slave  image.  Note  that  these  images  must  be  the  result  of  processing 
that  keeps  the  phases  intact,  so-called  phase-preserving  processing. 

We  denote  the  distance  from  the  source  to  the  target  rSJ  and  r{  and  r2  are  the  distances 
between  the  target  and  arrays  1  and  2,  respectively.  The  two-way  travel  times,  i.e. 
corresponding  to  the  path  from  the  source  to  the  target  and  back  to  array  1  or  2,  are 
denoted  by  ts\  and  /^.  These  are: 

*Sl=(rS+nk  {S2  =  (rS  +^2^.  (2-1) 

where  c  is  the  sound  speed.  In  practical  applications,  the  distance  between  the  source 
and  receivers  is  small  compared  to  the  distance  to  the  target.  The  time  axes  are  usually 
converted  to  the  average  distances  rS\  and  r $i\ 
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rS  [  ~ 


rS  +  rl 


rS  2  ~ 


rS  +  r2 


(2.2) 


such  that  the  travel  times  can  be  easily  converted  to  the  average  distances: 


ls\ 

2c 


rS  2  - 


lS2 

2  c 


(2.3) 


The  quantities  rs \  and  rsi  are  commonly  referred  to  as  the  range  axes  in  sonar  images. 
Image  display  with  respect  to  the  range  axis  (instead  of  the  two-way  traveltime) 
facilitates  the  interpretation  of  sonar  images. 


We  refer  to  the  image  constructed  with  the  data  recorded  on  array  l  as  the  master 
image,  the  image  corresponding  to  array  2  is  the  slave  image.  In  order  to  infer  the 
height  corresponding  to  the  image  pixel  (x,  rs i)  in  the  master  image,  we  have  to  find  the 
corresponding  image  point  in  the  slave  image.  If  we  succeed  to  determine  the  range 
difference  between  the  corresponding  pixels,  we  can  infer  the  height  of  the  pixel  w  ith  a 
geometrical  relationship.  It  should  be  explicitly  stated  that  in  the  scenario  under 
consideration,  path  differences  only  occur  at  the  detector  side.  As  a  consequence,  there 
is  a  factor  of  2  difference  between  the  one-way  range  difference  r2  -  rx  and  the  two-way 
average  range  difference  r &  -  rs \. 

Given  the  range  difference,  we  obtain  the  height  as  follows  (see  Figure  2.1).  First,  we 
estimate  the  look  angle  0  using 

cos(fl  +  C)=  r*~B  .  (2.4) 

V  2r{B 

In  this  expression,  <^is  the  angle  of  the  baseline  vector  with  the  vertical  axis  (tilt  angle). 
Then,  the  height  H  of  the  object  relative  to  a  reference  height  z-  0,  e.g.  the  seafloor  ,  is 
given  by  (see  also  Figure  2.1): 

"(nM  —  rx  cos  $  (2.5) 

Thus,  given  the  range  shift  between  corresponding  points  in  the  master  and  the  slave 
image,  we  can  determine  the  height. 

To  find  corresponding  points  in  the  master  and  slave  images,  first  a  co-registration  is 
carried  out.  This  step  uses  geometric  information  to  relate  events  detected  by  array  l  to 
those  detected  by  array  2.  After  this  step,  only  the  information  that  was  not  known  a 
priori  remains  to  be  retrieved.  This  task  is  easier  than  retrieving  all  information  from  the 
data,  i.e.  retrieving  already  known  information  plus  the  new,  yet  unknown  information. 
This  step  is  carried  out  in  step  two  by  estimation  of  the  instantaneous  phase  difference. 

Hence,  a  two-step  procedure  is  used  to  determine  the  range  shift: 

1  Image  co-registration; 

2  Estimating  instantaneous  phase  differences. 

In  the  image  co-registration  step,  a  range  shift  (which  may  vary  with  range)  is 
determined  using  information  about  the  geometry.  The  instantaneous  phase  information 
is  then  used  to  refine  the  resolution  of  the  height  information,  i.e.  to  obtain  local 
information  on  the  shape  of  objects  or  on  the  bathymetry. 
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2.2  Image  co-registration 


In  the  co-registration  step,  we  use  information  about  the  geometry  to  overlay  the  master 
and  slave  images.  This  is  achieved  by  computing  and  applying  a  correction  to  the  slave 
image  such  that  r\  =  r2  with  respect  to  the  reference  plane,  the  seafloor  (see  Figure  2.2). 
In  this  case  we  assume  that  the  seafloor  is  (more  or  less)  horizontal  and  that  we  know 
the  measurement  configuration,  i.e.  the  height  h  of  array  1  above  the  seafloor,  the 
baseline  B,  and  the  tilt  angle  are  known  (Figure  2.2).  In  that  case,  we  can  estimate  for 
each  range  rx  the  corresponding  range  r2  using 


(2.6) 


With  0c(r\ )  the  co-registration  look  angle  corresponding  to  the  direct  path  from  array  1 
to  the  seafloor.  The  co-registration  look  angle  0c{r\)  is  calculated  from: 


(2.7) 


Geometry-based  co-registration  thus  requires  information  about  the  acquisition 
geometry  and  a  rough  estimate  of  the  bathymetry.  When  this  information  is  not 
available  or  not  reliable,  data-driven  (correlation-based)  co-registration  methods  can  be 
considered.  These  have  been  thoroughly  investigated  by  [1]  and  [9]. 

When  r2  is  known,  we  can  determine  the  one-way  co-registration  range  shift  A rc: 


(2.8) 


•u source 


Z 


A  y 


Figure  2.2  Configuration  for  SAS  interferometry  (co-registration). 
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Implementation 


The  co-registration  procedure  comprises  two  steps: 

•  Compute  and  apply  the  co-registration  range  correction  A rc(n)  to  the  slave  image, 
i.e.  the  range  axis  (average  two-way  distance)  of  the  co-registered  slave  image 
Arc(r.) 


becomes: 


>52 


rS  2  — 


Since  this  range  shift  generally  differs  from  the  range  sampling,  the  master  and 
slave  images  sample  the  range  axis  (average  two-way  distance)  at  different  points. 
We  denote  the  difference  between  these  samples  by  Sr(rs i).  The  sampling  term 
yields  a  one-way  range  shift 


Arf(rSl)=2*KrSl) 


(2.9) 


We  need  to  take  this  term  into  account  since  we  did  not  resample  the  range  axis 
(i.e.  sampling  the  range  axis  at  identical  points  as  the  master  image)  after  the  co¬ 
registration  range  correction.  This  resampling  would  require  interpolation  starting 
from  a  non-uniformly  sampled  range  axis,  and  preserv  ing  the  instantaneous  phase  in 
this  interpolation  procedure  might  be  an  issue. 


•  The  total  one-way  range  difference  between  two  corresponding  pixels  contains  a 
contnbution  of  three  terms,  the  co-registration,  sampling  (both  determined  from  the 
co-registration),  and  the  range  difference  derived  from  the  instantaneous  phase: 

Ar,0,  =  Arc  +  Ars  +  Ar,  (2.10) 

In  the  next  section  we  discuss  the  instantaneous  phase  contribution  Ar,. 


2.3  Extracting  local  phase  information 


Image  co-registration  generally  has  a  low  resolution.  This  second  step  uses  local  phase 
information  extracted  from  the  master  and  slave  images  to  gain  resolution  in  the  height 
map.  This  step  is  thus  crucial  for  getting  information  on  the  shape  of  objects  or  on 
details  of  the  bathymetry.  The  analysis  of  complex  trace  attributes  enables  us  to  retrieve 
this  local  phase  information  [10].  This  can  be  measured  using  the  instantaneous  phase. 
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Figure  23  Complex-trace  analysis.  Real  (a)  and  quadrature  (b)  traces  for  a  portion  of  a  trace. 

The  envelope  is  shown  as  the  dashed  line  in  (a)  and  (b).  Instantaneous  phase  is  plotted  in  (c), 
reprinted  from  f  1 1  ] 


Figure  2.3  illustrates  complex  trace  analysis.  It  shows  a  real  and  the  corresponding 
quadrature  trace  (which  can  be  obtained  using  the  Hilbert  transform),  which  we  denote 
with  h r (/)  and  A /(f),  respectively.  The  complex  trace  h(t)  is  then 

h(t)  =  hr(t)+j.hi(t),  (2.11) 

where  j  denotes  the  complex  unit  vector.  The  envelope  is  defined  as  the  amplitude  of 
h(t ),  the  instantaneous  phase  is  given  by  the  phase  angle  of  h(t).  In  the  following,  we 
use  the  term  phase  information  for  the  instantaneous  phase.  For  band  limited  signals, 
the  phase,  or  the  phase  difference  A^can  be  related  to  the  travel  time  difference 
Ar=/-/0  according  to: 

A<f>  =  2/r/)Ar,  (2.12) 


Where /■  is  the  instantaneous  frequency  [10].  We  estimate  the  time  differences  and 
assign  these  to  the  centre  frequency  fc .  We  further  assume  that  the  centre  frequency 
fc~f .  Then,  the  instantaneous  phase  difference  corresponds  to  the  range  difference 


a  n  = 


cA(j> 

i^Tc- 


(2.13) 


As  illustrated  in  Figure  2.4,  instantaneous  phase  gives  only  accurate  time  or  range 
differences  for  narrowband  data  (and  for  small  time  shifts).  Figure  2.4  shows  that  for  a 
broadband  LFM  pulse  (fc  =  150  kHz,  b  -  60  kHz),  the  error  is  30%  in  the  estimated  time 
differences.  Therefore,  it  is  worthwhile  to  consider  data  analysis  in  different  frequency 
bands.  As  an  alternative,  correlation  methods  can  be  used.  The  discontinuous  jumps  in 
the  curves  are  a  result  of  cycle  skips.  Figure  2.4  also  illustrates  the  importance  of  the 
image  co-registration.  Instantaneous  phase  estimates  are  more  accurate  for  small  time 
differences.  Finally,  the  right-hand  graph  shows  that  this  type  of  error  is  not  an  issue  for 
radar  interferometry,  because  the  relative  bandwidth  is  usually  much  smaller. 
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Figure  2.4  Difference  between  the  ‘true  time  shift’  At  and  the  time  shift  derived  from  the  instantaneous 
phase  as  a  function  of  the  time  shift  for  matched  filtered  LFM  pulses  with  150  kHz  centre 
frequency  and  different  bandwidths.  The  right  graph  shows  this  sensitivity  analysis  for  a  radar 
scenario  with  a  15  GHz  centre  frequency 


2.4 


Implementation 


We  denote  the  co-registered  images  by  U\(m,n)  and  Ui{m,n)y  where  the  indices 
m  =  1,  . ..,  rtr  and  n  =  1,  . nx  denote  the  samples  of  the  range  and  azimuth  axis, 
respectively.  Both  U\(m,n)  and  Uiimji)  are  complex.  Interferometry  uses  the 
differences  in  the  instantaneous  phase  to  obtain  information  on  the  height  of  the  imaged 
area.  The  procedure  consists  of  the  following  steps: 


Measure  phase  difference  and  construct  the  interferogram; 
We  first  define  the  correlation  map  C(m,ri): 


k  / 


C(w,h)= 


k  l 


X  + P<n  +  q)j*1(m  +  p,n  +  q) 

p=-k  g=-l _ 

v/v  *  ,  v/2 


(2.14) 


X  YP'(m+ P'n+A~  X  XI \ui(m+ P'n+<}\ 

\p~-k  q=-l  )  \P=~k  q--l 


The  correlation  map  is  computed  averaging  over  patches  w  ith  a  size  of  2A+1  times 
2/+1  elements  The  amplitude  of  C(m,fl)  is  considered  to  be  an  important  measure 
of  data  quality  [13],  and  is  often  used  for  creating  a  mask.  Nevertheless,  this 
measure  is  rather  sensitive  to  phase  differences  between  the  master  and  slave 
images.  Therefore,  we  use  the  correlation  map  of  the  instantaneous  amplitudes  as  a 
measure  for  data  quality  instead. 


The  instantaneous  phase  of  C(m,/f)  yields  the  phase  difference  between  the  master 
and  slave  images  corresponding  to  image  pixel  {m,n).  The  2-D  moving  average 
(summation  of  neighbouring  pixels)  is  applied  to  improve  the  signal-to-noise-ratio 
(SNR)  of  the  phase  estimate.  Note  that  filtering  of  the  correlation  map  is  preferred 
over  filtering  of  the  phase  values  only  [13] 


•  Apply  phase  unwrapping; 

Phase  unwrapping  is  required  since  the  instantaneous  phase  is  bounded  between  -i r 
and  7r  (see  Figure  2.3).  As  a  result,  ambiguity  exists  in  the  phase  values,  since  a 
multiple  factor  of  2 n  is  left  undetermined.  Many  algorithms  exist  to  unwrap  the 
phase.  These  algorithms  are  based  on  minimizing  discontinuities  in  the  data. 
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minimizing  the  length  of  branch  cuts,  and  least  squares  algorithms.  For  a  review  of 
phase  unwrapping  methods,  we  refer  to  [13].  In  our  current  implementation,  we  use 
the  classical  branch-cut  (Goldstein)  algorithm  [12],  [13].  2-D  phase  unwrapping  is 
an  active  field  of  research.  For  further  development  it  would  be  worthwhile  to 
investigate  also  so-called  "net^'ork  approaches'  and  ‘ synthesis  algorithms'  [14], 
[15],  [16]. 

•  Estimate  the  height; 

First  estimate  the  range  difference  A rtot  using  Equation  2.10,  then  determine  r2\ 

r2=ry+  Arlol .  (2.15) 

The  value  of  r2  is  now'  used  to  compute  the  angle  according  to  Equation  2.4,  after 
which  the  height  can  be  obtained  using  Equation  2.5. 
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3  Interferometric  SAS  processing:  simulated  examples 


We  proceed  by  illustrating  SAS  interferometry  with  simulated  examples. 

Before  discussing  the  examples,  we  first  review  the  essential  features  for  the  simulation 
and  processing  [17],  and  analyse  whether  the  phase  information  is  preserved. 


3.1  SIMONA  SAS  Simulator 

In  this  section  we  briefly  review  the  main  functionality  of  the  SIMONA  simulator. 

Let  xr  and  xs  denote  the  receiver  and  source  positions  and  s(/)  the  source  signal, 
respectively.  The  simulator  models  the  wavefield  U(t ,  xr ;  x*)  scattered  by  a  set  of  N 
point  scatterers.  It  is  given  by: 

N 

U{xr;xs,t)  =  -  T^x^x,)) .  (3.1) 

1=1 

For  each  point  scatterer,  the  response  is  obtained  by  convolving  (*)  the  source  signature 
with  the  target  response.  It  is  characterized  by  a  time  delay  r(xr;  x^)  and  an  amplitude 
response  A(xr\  x^).  The  time  delay  is  the  propagation  time  from  the  source  to  the  target 
and  back  to  the  receiver.  The  amplitude  term  includes  both  target  strength  and 
propagation  loss. 

Equation  3. 1  is  evaluated  at  discrete  samples  t  =  kAt : 

N 

U{\ r ;  x 5 , kAt )  =  s(kAt) *  ^  At  (x r ;  x , ]S(kAt  -  r,  (x r ; x , ))  (3.2) 

/=1 

Since  the  time  delays  rz(x,;  xs)  do  not  necessarily  coincide  with  the  discrete  time 
samples  kAt ,  the  wavefield  obtained  with  Equation  3.2  does  not  guarantee  that  the  phase 
information  is  correctly  represented  in  U(xr ;  xs,  kAt).  Recall  that  phase  errors  can  be 
substantial  when  the  wavefield  is  sampled  according  to  the  Shannon  sampling  criterion 
(at  least  two  samples  are  required  per  minimum  wavelength). 

In  order  to  preserve  the  phase  information,  we  evaluate  the  response  of  the  scattered 
field  in  the  frequency  domain.  Defining  the  Fourier  transform  (FT)  as 

OD 

/(«>)=  |/(?)exP[-  (3.3) 

—00 


with  inverse  transform  (I FT) 

1  °° 

fit)  =  —  (f{(o)exv[ia>t)da> , 


Equation  3.2  can  be  written  as: 


£/(xr;x5,£A/)  =  IFT 


N 

s{kA(o)^  Ai{xr;xs  )exp[-  ikAeorl  (x  r ;  x  s )] 


i=i 


(3.4) 


(3.5) 
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In  this  equation,  the  exact  values  for  the  time  delays  are  used.  As  a  result,  phase 
information  is  correctly  captured  in  U(xr ;  x*  kAt).  The  implementation  for  target 
strength  is  identical  to  the  existing  SAS  Simulator  (SIMONA  [17]). 

At  this  point  it  has  to  be  stated  that  the  existing  SAS  simulator  (SIMONA  [17])  did  not 
capture  the  phase  information  correctly.  Phase  errors  (on  the  order  of  A732)  were 
introduced  in  a  procedure  to  optimise  SIMONA  for  speed.  It  has  been  verified  that  these 
phase  errors  do  not  influence  the  SAS  imaging.  The  reason  is  that  these  images  display 
the  instantaneous  amplitude.  As  shown  in  Figure  2.3,  the  instantaneous  amplitude  varies 
smoothly  in  comparison  to  the  instantaneous  phase,  and  is  therefore  not  sensitive  to 
small  phase  errors.  Interferometry,  on  the  other  hand,  is  based  on  the  interpretation  of 
subtle  phase  differences,  and  these  phase  differences  may  therefore  be  significant  to 
interferometry  [3].  The  phase  errors  introduced  (rounding  errors)  differ  for  each  point 
scatterer  and  also  for  each  array  element. 

For  the  development  of  interferometry,  we  used  the  frequency-domain  implementation 
described  in  this  section.  This  procedure  guarantees  that  the  phase  information  is 
correctly  captured  in  the  simulated  data. 

Software 

The  SAS  simulator  (Sourceforge  “)  Source  Code  ->  SAS  Simulator )  has  the  following 
structure: 

•  create_input\main_input_sassimulator.m:  this  code  enables  the  user  to  specify  the 
input  to  the  SAS  simulation  and  saves  these  settings  in  a  parameter  file  In  the 
examples  presented  here,  different  targets  have  been  used,  and  these  are  referred  to 
with  a  character  string,  e.g.  id  =  'testlO';  create  Jnput\define_target_testl0.m. 

•  simulator\main_simulator.m  computes  the  pingfiles  for  the  specified  input  and 
saves  these  in  a  directory  outputLtimetag  (yyyymmddhhmm). 

In  the  following,  we  refer  to  the  simulation  results  using  the  character  string  Id. 

3.2  SAS  Processing:  Imaging  using  Stolt  migration 

Imaging  or  migration  is  the  procedure  that  converts  the  zero-offset  (phase  centre) 
section  of  the  measured  data  into  an  image  of  the  medium.  The  following  procedure 
assumes  that  the  zero-offset  section  is  a  uniform  linear  array,  and  we  denote  the 
midpoint  coordinate  x.  This  coordinate  is  the  average  of  the  source  and  receive  element 
positions. 

As  in  the  papers  describing  the  Stolt  migration  [17],  [18],  we  distinguish  the  recorded 
data  from  the  migrated  data  introducing  the  arguments  U(r  =  0,  x ,  /)  to  denote  the 
recorded  data  and  U (r,  xy  t  =  0)  for  the  imaged  data,  respectively.  The  recorded  data  are 
acquired  at  location  x  in  the  form  of  time  series,  and  migration  converts  these  time 
series  into  an  image  of  the  medium  itself.  This  is  represented  by  the  scattering  strength 
as  a  function  of  location  (rrt),  where  r  indicates  the  range.  The  migration  procedure 
converts  hyperbolas  to  point  scatterers  corresponding  to  the  correct  image  positions. 

For  underwater  imaging,  the  migration  can  be  efficiently  implemented  using  Stolt 
migration  [17],  [18],  [19].  This  method  assumes  that  the  sound  velocity  in  the  medium 
is  constant,  which  may  lead  to  some  defocusing  in  case  of  propagation  through 
sediment  (for  buried  objects).  For  now,  we  assume  that  this  defocusing  effect  is 
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negligible,  since  we  mainly  consider  shallowly  buried  objects  and  steep  incidence 
angles.  Additionally,  the  implementation  requires  a  uniform  linear  data  array. 
Stolt  migration  evaluates  U ( r ,  x,  t  =  0)  in  the  frequency-wavenumber  domain: 


U(r,x,6)  =  — —  \dco  \dkx  U(0,kx,co)exp 

2n  J  J 


4(0 z  ,2  , 

— — kx  r  +  kxx 


(3.6) 


where  U (0,  kx,  co)  is  obtained  by  taking  the  Fourier  transform  of  the  measured  data 
U(0,x,t), 


OO  cojc 

U(0,kx,co)  =  —  [ dt  [cir  (/(0,x,/)exp[-  i(cot  +  A:x.r)] 
2 k  J  J, 

-GO  -cojc 


(3.7) 


Stolt  migration  is  based  on  the  change  of  variables  from  the  angular  frequency  co  to  the 
wavenumber  kr  according  to  the  formula: 


In  order  to  compute  Equation  3.6,  we  apply  the  coordinate  transform  and  change  the 
integration  variable: 


00  cojc 

U(r,x, °)  =  —  } dkr  | dkx 

-oo  -cojc 


t/(/:r(dio),/:JC,0)exp[/(/:r  r  +  kxx)\. 


(3.9) 


The  evaluation  of  the  Fourier  transform  (with  respect  to  kr)  requires  that  the 
wavenumber  is  discretised  using  equal  spacing.  This  implies  that  we  need  U( 0,  kx>  co)  at 
non-equally  spaced  locations  which  are  obtained  from  Equation  3.8.  These  values  are 
obtained  by  interpolation  (spline). 


Implementation 

•  2-D  Fourier  transform  of  recorded  data. 

•  Apply  phase  shift  to  account  for  time  window  to :  multiply  U( 0,  kx ,  co)  with 
exp[/a>/0]; 

Shifting  property  :  FT (/(())  ++f((o)\  FT (f(t  +  t0))  <->  f(a>)  exp[/£U/0] 

Input  to  Fourier  transform:  f(t -to),  which  corresponds  to  multiplying  with 
exp[-/(y/0]  in  the  frequency  domain. 

•  Multiply  now  with  exp [ikrrQ  ]  ,  where  kr  is  given  in  Equation  3.8.  This  step  has  to 
be  taken  before  the  interpolation,  and  define  the  corresponding  range  axis. 

•  Prepare  for  inverse  Fourier  transform:  variable  substitution  and  computation  of 
amplitude  term  (Jacobian). 

•  Interpolation  to  get  U{kr{co\kxfi]  uniformly  sampled.  Approach:  Define  a  linear 
array  for  kr  and  compute  the  corresponding  frequency  values  using  Equation  3.8. 
Repeat  this  step  for  each  value  of  the  horizontal  wavenumber.  Use  spline 
interpolation  to  obtain  U{kr  (&>),  kx  ,0) . 

•  2-D  Inverse  Fourier  transform. 
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Software 

The  SAS  processor  (Sourceforge  ■)  Source  Code  ->  SAS_Processor)  is  organised  as 
follows: 

•  The  folder  dataQC  contains  scripts  for  data  visualisation  and  quality  control. 

•  The  folder  datalO  contains  scripts  for  data  input-output. 

•  The  preprocessing  script  (i.e.  create_sasinput_DERAGESMA_raildata_051003.m) 
constructs  the  synthetic  array  from  the  individual  ping  files.  The  synthetic  array  is  a 
uniform  linear  array.  Motion  compensation  is  an  essential  step  in  the  preprocessing 
procedure. 

•  The  folder  imaging  contains  the  migration/imaging  codes 
(i.e.  create_sasoutput_simulations.m). 

3.3  Example  I:  tilted  plane7 

The  first  simulated  example  concerns  a  tilted  plane,  as  illustrated  in  a  2D  and  a  3D  view 
in  Figure  3.1.  The  target  is  built  up  from  individual  point  scatterers  that  are  separated  by 
X/5  =  3.0  mm.  The  target  roughness  is  also  3.0  mm.  The  simulation  parameters  are 
listed  in  Table  3.1. 


41J 


Figure  3  1  Top  view  (left)  and  3-D  view  (nght)  of  the  tilted  plane 


Table  3. 1  Simulation  parameters  for  the  tilted  plane  example. 


100  kHz 
40  kHz 


Pulse  centre  frequency  fc 

Pulse  bandwidth  b 

Pulse  type 

Pulse  duration 

Length  of  synthetic  aperture 

Number  of  pings 

Number  of  elements/ping 

Separation  of  elements 

Displacement  between  pings 

Platform  height 

Platform  y-distance  to  target 

Number  of  receive  arrays 

Baseline  vector  B  (vertical  separation) 

Source  position  _ 


Linear  frequency  modulated  pulse  (LFM) 
6  4  ms 

12  m  (from  0  0  to  12.0  m) 

200 

8 

0.015  m 
0.060  m 
10.0  m 
40  0  m 
2 

(0,  0t  0  20)  m 

Centre  of  array  1 _ 


i 
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Having  simulated  the  data,  we  apply  Stolt  migration  to  convert  the  SAS  data  to  the 
image  of  the  reflectivity  of  the  tilted  plane  (Figure  3.2).  Note  that  we  did  not  include  a 
background  in  the  simulation. 
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Figure  3,2  Simulated  synthetic  aperture  sonar  data  (left)  and  migrated  data  (nght)  of  the  tilted  plane 
illustrated  in  Figure  3.1,  No  background  noise  was  simulated. 


Subsequently,  the  height  map  is  derived  using  interferometric  processing: 

•  images  are  constructed  using  the  data  measured  on  both  arrays; 

•  geometry-based  co-registration  is  applied  with  respect  to  2  =  0; 

•  instantaneous  phase  differences  are  estimated; 

•  2-D  phase  unwrapping; 

•  height  estimation. 

Figure  3.3  shows  details  of  the  SAS  images  derived  from  both  arrays.  The  range  shift 
induced  by  the  baseline  is  approximately  0.02  m.  The  geometry-based  co-registration 
compensates  for  this  range  shift.  Similar  to  the  configurations  shown  in  Figure  2.1  and 
Figure  2.2,  array  1  is  the  low  er,  array  2  the  upper  of  the  tw  o.  Therefore,  rz  >  r  1,  w  hich 
explains  the  0,02  m  shift  in  range  between  the  left  and  the  centre  panel  of  Figure  3.3. 


Figure  3.3  Detail  of  SAS  images  constructed  from  the  lower  array  1  (left  panel)  and  the  upper  array  2 
(centre  panel)  The  geometrical  range  displacement  between  the  images  is  approximately 
0  02  m  The  nght  panel  shows  the  co-registered  slave  image. 


Figure  3.4  displays  the  height  estimated  using  the  interferometric  processing.  For  image 
pixels  with  a  low  scattering  strength  (-40  dB  relative  to  the  maximum  value),  we  set  the 
height  to  2  =  0  (the  phase  information  is  then  unreliable).  A  comparison  of  the 
interferometric  height  estimate  to  the  shape  of  the  3-D  target  (Figure  3.1)  shows  a  good 
agreement.  This  is  confirmed  by  the  3-D  reflectivity  image  (Figure  3.5).  This  image 
displays  the  reflectivity  as  a  function  of  azimuth,  range,  and  height. 

The  noise  on  the  height  estimates  is  a  concern.  Figure  3.5  shows  that  the  peaks  and 
valleys  in  the  height  estimates  generally  have  a  low  reflectivity  strength.  This  indicates 
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that  sidelobes  from  strong  scatterers  (highlights)  might  influence  the  phase  information. 
This  results  in  unreliable  phase  information  since  it  is  not  related  to  the  corresponding 
image  position  in  such  a  situation.  Similar  reasoning  applies  to  the  height  estimates  for 
locations  outside  the  target. 


A™“*w  wwi'H 

Figure  3.4  Interferometric  height  estimate  (left)  and  target  (right). 


Figure  3,5  3-D  view  of  reflectivity  (left)  and  3-D  target  (right).  Note  that  the  peaks  and  valleys  in  the 

reflectivity  image  have  a  low  reflectivity  strength. 


3.4  Example  II:  top  part  of  a  sphere2 

In  the  second  simulated  example  we  consider  the  top  part  of  a  sphere  as  a  target. 
The  centre  of  the  hemispherical  target  is  located  at  (x,>\  z)  =  (6,  40,  0)  m,  see  also 
Figure  3.9,  right  panel.  The  radius  of  the  full  sphere  is  0.35  m.  The  other  simulation 
parameters  are  identical  to  the  previous  example  (Table  3.1). 


Figure  3  6  Simulated  synthetic  aperture  sonar  input  data  (left)  and  migrated  data  (right)  for  array  1 


2 


Simulations,  processing  and  interferometry-  Software  id  ‘test7_2 0080 128’. 
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Figure  3.6  shows  the  SAS  input  data  and  the  corresponding  image  for  array  1 . 
Characteristic  in  images  of  a  sphere  is  the  strong  highlight  induced  by  normally  incident 
waves;  the  target  back-scattering  strength  is  large  for  waves  arriving  at  normal 
incidence.  Diffuse  scattering  is  determined  by  target  roughness  at  different  angles  of 
incidence,  and  is  therefore  significantly  smaller. 

Figure  3.7  shows  details  of  the  SAS  images  derived  from  both  arrays.  The  range  shift 
induced  by  the  baseline  is,  similar  to  the  previous  example,  approximately  0.02  m. 

The  geometry-based  co-registration  (with  reference  level  z  -  0)  compensates  for  this 
range  shift.  Note  that  sidelobes  occur  in  the  range  direction  due  to  the  high  amplitude 
signal. 


Figure  3.7  Detail  of  SAS  images  constructed  from  the  lower  array  1  (left  panel)  and  the  upper  array  2 
(centre  panel).  The  geometrical  range  displacement  between  the  images  is  approximately 
0.02  m  The  right  panel  shows  the  co-registered  slave  image 


Figure  3.8  shows  the  height  inferred  from  interferometry.  For  image  pixels  with  a  low 
scattering  strength  (-40  dB  relative  to  the  maximum  value),  we  set  the  height  to  zero. 


Figure  3.8  Interferometric  height  estimate  for  the  hemispherical  target.  Note  that  the  radius  is  0  35  m. 


Figure  3.9  3-D  view  of  reflectivity  (left)  and,  for  comparison,  the  3-D  target  (right). 


TNO  report  |  TNO-DV  2008  A176 


22/44 


We  observe  in  Figure  3.9  that  the  height  increases  towards  the  centre  of  the  sphere. 

The  shape  of  the  sphere  is  quite  well  reconstructed.  Nevertheless,  we  do  not  fully 
recover  the  shape  of  the  sphere.  Two  mechanisms  are  identified  which  may  explain  this: 

•  The  scattering  strength  of  the  highlight  is  dominant.  As  a  result,  sidelobes  may 
influence  the  phase  esttmate  at  other  points  of  the  sphere. 

•  There  is  height  ambiguity  in  range  direction,  as  illustrated  in  Figure  3.10. 

The  first  response  of  a  spherical  or  cylindrical  target  is  recorded  for  a  wave  at 
normal  incidence.  When  the  range  increases  with  8r,  signal  is  backscattered  from 
two  different  heights  and  is  thus  recorded  simultaneously  on  array  1 .  Only  the  path 
lengths  to  array  2  differ  slightly.  As  a  result  of  this  range  ambiguity,  the  time  (or 
phase)  differences  to  the  individual  scatterer  positions  cannot  be  resolved. 

This  mechanism  is  illustrated  in  Figure  3,11.  It  illustrates  that  we  only  obtain  an 
estimate  of  the  phase  difference  corresponding  to  the  average  height. 


Figure  3.10  Height  ambiguity  with  range;  the  first  response  of  a  spherical  or  cylindrical  target  is  recorded 
for  a  wave  at  normal  incidence.  When  the  range  increases  with  signal  is  backscattered  at 
two  different  heights  and  is  thus  recorded  simultaneously  on  array  1.  The  path  lengths  of  these 
two  scatterer  positions  to  array  2  differ  slightly. 


Figure  3.11  Conceptual  illustration  of  effect  of  range  ambiguity  on  recorded  waveforms  and  on  the 
estimated  time  difference  The  waveform  recorded  on  array  1  at  range  r  and  the 
corresponding  waveform  on  array  2  (left  panel)  The  additional  path  lengths  of  the  individual 
scatterers  to  array  2  are  0.05  and  0.20  s,  respectively.  The  cross-correlation  of  the  signals 
received  at  array  1  and  2  shows  that  only  the  average  of  these  path  length  differences  can  be 
resolved  (right  panel). 
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4  SAS  interferometry  applied  to  DERA/GESMA  rail  data 


4.1  Introduction  to  the  DERA/GESMA  rail  datasets 

In  1999,  rail  based  experiments  were  conducted  by  DERA  (Defence  Evaluation  and 
Research  Agency  -  UK)  and  GESMA  (Groupe  d’Etudes  Sous-Marines  de  FAtlantique  - 
FR)  in  the  bay  of  Brest  at  Lanveoc  (FR).  The  goal  of  this  campaign  was  to  gather  data 
in  a  semi-controlled  environment  to  further  improve  different  types  of  algorithms  such 
as  SAS  and  InSAS. 

The  sonar  system  used  is  mounted  on  a  rail  (see  Figure  4.1)  and  is  composed  of  tw  o 
receiving  arrays  of  32  elements  each,  vertically  separated  by  21.7  cm.  A  source  is 
placed  under  the  two  receiving  arrays,  that  emits  LFM  pulses  in  the  frequency  band  of 
120  kHz  to  180  kHz.  Detailed  parameters  of  the  acquisition  are  presented  in  Table  4.1. 


Table  4  1  Acquisition  parameters. 


Parameters  of  the  system 

Value 

Centre  frequency  [kHz] 

150 

Sampling  frequency  [kHz] 

560 

Bandwidth  [kHz] 

60 

Pulse  length  [s] 

0,004 

Pulse  repetition  time  [s] 

3.34 

Platform  speed  [m/s] 

0  01 

Space  between  elements  [m] 

0  00834 

Displacement  between  pings  [m] 

0.0334 

Total  number  of  elements 

64 

Number  of  elements  per  array 

32 

Vertical  array  separation  [m] 

0.217 

Total  number  of  pings 

320 

Sound  velocity  [m/s] 

1495 

Number  of  samples 

16384 

On  the  test  field  across  the  test  rail,  several  mine-like  objects  were  placed  on  the  sea 
floor  at  four  different  ranges.  A  sketch  and  corresponding  sidescan  survey  of  the  test 
area  are  found  in  Figure  4.2  and  Figure  4.3,  respectively. 
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Figure  4.2  Rail  experiments  (GLSMA-DLRA,  1999)  The  sonar  system  is  mounted  along  the  quay  and 
faces  the  measurement  field.  The  target  field  is  divided  in  25  m  zones  where  different  targets 
are  deployed.  At  15  m,  different  targets  are  buried  a  sphere,  and  two  mine-like  targets:  a 
rockan  and  a  manta.  At  20  m  a  mine-like  sphere  is  also  buned  At  25  m  a  mine-like  sphere  is 
placed  on  the  seabed  Three  mine-like  targets  are  placed  at  50  m  range  a  manta,  a  rockan  and  a 
sphere 

At  75  m  two  more  mine-like  targets:  a  cylinder  and  a  sphere  and  at  100  m  a  mkl2  mine  like 
cylinder.  The  corresponding  datasets  are  mentioned  at  the  right  side  of  the  field 
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Figure  4  3  Side  scan  survey  of  the  test  area  on  the  left  we  can  distinguish  the  rail  where  the  sonar  was. 

On  the  right  side  of  the  rail,  along  the  y-axis,  the  targets  at  three  different  ranges  (25,  50,  75  m) 
can  be  distinguished. 
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Table  4,2  Experiment  descriptions. 


Experiment  Number 

Description 

051002 

Spherical  mine-like  target  at  25  m 

051003 

011008 

Manta,  spherical  and  rockan  mine-like  targets  at  50  m 

011004 

Spherical  and  cylindrical  mine-like  targets  at  75  m 

The  different  tested  targets 

are  mine-like  objects. 

Table  4  3  Targets. 

Target 

Description 

Sphere 

1  m  diameter 

Cylinder 

55  cm  diameter  and  90  cm  long 

Despite  the  parameters  given  in  Table  4.1,  a  few  parameters  are  still  unknown  and  the 
most  important  ones  are  the  platform  height  above  the  seafloor,  which  is  one  of  the 
required  parameter  for  InSAS  ( h  in  Equation  2.5)  and  the  range  vector  corresponding  to 
the  data  acquisition  (i.e.  a  time  reference).  These  parameters  had  to  be  deduced  and 
estimated  from  the  acoustical  data  itself,  because  they  were  not  documented.  In  order  to 
determine  the  range  bracket,  a  SAS  migration  for  various  range  vectors  was  performed 
trying  to  optimise  the  contrast  on  the  spherical  target.  An  estimate  of  the  height  of  the 
platform  is  inferred  from  the  shadow  length.  The  parameters  for  experiments  0 1 1 004 
and  051003  are  summarised  in  Table  4.4  and  Table  4.5,  respectively. 

4.2  Description  of  processing  sequences 

4.2. 1  Quality  control 

Before  any  SAS  processing  is  applied  to  the  raw  data,  it  has  to  undergo  a  quality  control 
sequence  in  several  steps  in  order  to  optimize  the  final  SAS  picture.  For  further 
reference,  channels  1  to  32  correspond  to  the  lower  array  (later  referred  to  as  array  1) 
and  channels  33  to  64  correspond  to  the  upper  array  (later  referred  to  as  array  2). 


(a)  (b) 


Figure  4  4  Ping  1 60.  Raw  data  (a)  and  the  same  dala  after  amplitude  equalization/inlerpolation  (b). 

A  limited  range  bracket  from  one  ping  of  the  01 1004  dataset  is  presented  in  Figure  4.4. 
The  left  panel  (a),  shows  that  defective  channels  occur,  containing  invalid  data. 

This  may  be  due  to  bad  connections  and  short-circuits  in  the  array  electronics,  causing 
duff  channels  (blue  columns  containing  zeroes)  and  -less  visible-  identical 
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neighbouring  channels,  respectively.  These  defective  channels  are  restored  by 
interpolation  of  the  neighbouring  channels,  resulting  in  Figure  4,4  (b). 

Additional  to  channel  restoration,  amplitude  equalisation  is  performed.  Channel 
amplitudes  may  vary  slightly  due  to  differences  in  element  sensitivity,  gain  factors  and 
other  causes,  as  illustrated  in  Figure  4.5. 

In  order  to  correct  these  small  effects,  amplitude  equalization  is  performed  on  each  ping 
file:  for  each  ping  file,  the  total  energy  Et  along  a  channel  i  is  computed: 

N 

Ex  -  (/)  ,  where  N  is  the  number  of  samples.  Then  each  trace  i  is  normalized 

M 


Figure  4.5  Average  energy  per  channel  in  red  and  standard  deviation  in  black.  This  figure  is  an  average 
across  all  pings  for  dataset  01 1004. 


4.2.2  Synthetic  aperture  formation 

In  order  to  allow  coherent  SAS  migration,  an  equivalent  array  of  midpoints  is  used 
instead  of  the  physical  receive  element  positions.  These  monostatic  midpoints  are 
located  in  the  middle  between  each  element  and  the  position  of  the  acoustic  centre  of 
the  source  (see  Figure  4.6).  Hence  the  effective  array  length  is  only  half  the  physical 
array  length. 
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Figure  4.6  Phase  centre  approximation:  each  transmitter-receiver  pair  is  assumed  to  be  a  monostatic 
equivalent  located  midway  between  the  two. 

The  integration  of  successive  pings  is  then  done  for  common  midpoints,  i.e.  midpoints 
of  successive  pings  that  are  at  the  same  location  (Figure  4.7). 


Azimuth  [m] 


Figure  4.7  Four  consecutive  pings  (160  to  163):  for  clarity,  the  arrays  are  shifted  in  range  in  this  figure. 

The  circles  represent  the  actual  element  locations  and  the  pluses  (+)  are  the  corresponding 
midpoint  locations.  The  common  midpoints,  i  e.  with  the  same  position,  are  circled  in  red  for 
one  azimuth  position  The  signals  from  these  four  overlapping  midpoints  are  summed. 

In  order  to  check  if  overlapping  midpoints  are  indeed  in  phase,  it  is  possible  to  check 
visually  if  the  corresponding  traces  for  one  position  are  in  phase  or  not. 
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Figure  4.8  Common  midpoints  stacks  zoomed  in.  The  first  four  curves  correspond  to  the  response  on  the 
following  midpoints:  l)ping  147,  element  31,2)  ping  148,  element  23,  3)  ping  149,  element 
1 5,  4)  ping  150,  element  7.  The  right-curve  is  the  normalised  summation  of  the  four  common 
midpoints.  In  theory  the  five  curves  should  be  exactly  aligned  in  range  and  with  the  same 
amplitude,  which  is  not  the  case;  there  are  indeed  some  (small)  shifts,  and  as  a  result  the 
summed  curve  has  a  lower  amplitude 

The  misalignment  visible  in  Figure  4.8  is  caused  by  an  offset  in  the  position  of  the 
receive  arrays  themselves.  Despite  being  on  a  rail,  it  was  found  that  the  sonar  system  is 
not  perfectly  aligned  with  the  rail  axis,  leading  to  a  crabbing  angle  (Figure  4.9). 

Other  motions  (especially  sway  -sideways  displacement-  and  yaw  -rotation  around  the 
vertical  axis-)  are  negligible  due  to  the  nature  of  the  platform.  This  position  error  must 
be  corrected  to  allow  a  coherent  summation  of  the  different  traces. 


Antenna 

Antenna 

Antenna 

Antenna 

Pingn 

Pingn+1 

Pingn+2 

Ping  n+? 

Direction  of  motion 


Figure  4.9  Illustration  of  magnified  crabbing  angle,  constant  angle  between  the  axis  of  the  antenna  and  the 
direction  of  motion 
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4.2.3  Motion  compensation 

A  rail  is  an  extremely  stable  platform  in  terms  of  motion,  especially  compared  to  a  ship. 
Still,  for  SAS  imaging,  a  motion  error  exceeding  A/10,  where  A  corresponds  in  this  case 
to  the  highest  frequency  transmitted,  could  cause  defocusing  of  the  output  image. 

This  corresponds  to  an  error  greater  than  0.8  mm  in  our  case,  which  is  very  small. 
Moreover,  SAS  interferometry  is  even  more  sensitive  to  motion  errors  than  SAS 
imaging.  The  rail  can  therefore  not  be  considered  as  a  linear  track  and  its  accurate 
motion  has  to  be  derived  from  the  acoustic  signals  themselves.  This  approach  is 
generally  referred  to  as  autofocusing.  As  mentioned  previously,  a  crabbing  angle  was 
found  in  the  system,  which  has  to  be  corrected.  All  other  motions  (especially  yaw  and 
sway,  which  are  most  disturbing  for  SAS  imagery)  were  found  negligible. 

A  well  known  and  used  technique  for  this  purpose  is  the  displaced  phase  centers  (DPC) 
method  [17],  [20].  This  method  uses  the  previously  mentioned  phase  centers  and  is 
based  on  correlating  two  consecutive  pings:  A  matrix  of  absolute  cross-correlation 
values  can  be  created  for  tw  o  consecutive  pings,  consisting  of  all  possible  combinations 
of  array  elements.  The  highest  values  in  the  matrix  indicate  the  displacement  between 
those  pings,  expressed  in  number  of  array  elements.  The  example  in  Figure  4. 10  shows 
a  displacement  of  8  elements  along  the  array  axis  (along  track). 


Figure  4. 1 0  Displaced  phase  centre  method  for  motion  estimation*  on  the  X-axis  there  is  ping  n  and  on  the 
Y-axis  there  is  ping  (n+1)  The  represented  values  correspond  to  the  correlation  coefficients. 

The  current  DPCA  algorithm  cannot  distinguish  between  a  fixed  crabbing  angle  and 
constant  sway  as  illustrated  in  Figure  4. 1 1 .  It  illustrates  the  correspondence  between  an 
angle  corresponding  to  a  constant  sway  drift  and  a  constant  crabbing  angle: 

®sway  -  @  crabbing  •  ^ this  ls  interpreted  as  a  sideways  displacement,  this  will  result  in 
increasing  lateral  displacements  along  the  path,  as  shown  in  Figure  4.12. 

For  our  purpose,  i.e.  estimating  the  crabbing  angle,  the  sway  displacement  had  to  be 
derived  first  before  being  converted  to  a  crabbing  value.  First,  the  along  track 
displacement  is  estimated  from  the  displacement  between  the  main  diagonal  and  the 
diagonal  with  the  highest  values  (see  Figure  4.10).  This  gives  the  corresponding  phase 
centers  for  two  consecutive  pings.  The  time  lag  of  the  cross-correlation  peak  then  yields 
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the  displacement  in  the  range  direction  {cross-track).  The  average  value  of  all 
overlapping  elements  yields  the  sway  estimate. 


Configuration ;  SWAY 
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Figure  4  1 1  Top  figure  is  the  configuration  of  a  sway  motion,  the  DPC  algorithm  assumes  that  we  have  this 
situation  and  between  each  ping  it  computes  the  sway  Ay.  But  in  case  of  rail  measurements  it  is 
not  possible  to  have  such  a  configuration.  However,  by  rotating  the  top  figure,  we  obtain  the 
bottom  one  that  represents  the  crabbing  angle  configuration 


Conf  Iguralionj.  CRABBING 


The  DPC  algorithm  is  thus  applied  to  the  total  considered  dataset  and  the  trajectory'  of 
each  channel  is  estimated  (in  blue  for  only  one  element  in  Figure  4.12).  The  green  line 
is  the  linear  regression  corresponding  to  the  sway  displacement  and,  as  just  explained, 
this  gives  us  a  direct  estimate  of  the  crabbing  angle. 


Correction  of  the  trajectory-  hydrophone  16 


Figure  4. 12  Estimation  of  the  crabbing  angle,  only  element  16  is  represented  so  that  the  draw  ing  is  clearer. 
In  blue  the  lateral  displacement  estimated  w  ith  the  DPC  algorithm  and  in  green  the 
corresponding  linear  regression.  In  red,  the  projection  of  the  blue  curve  on  the  along  track,  i.e. 
the  real  movement  of  the  sonar  and  the  black  line  is  the  rail  axis. 

Once  the  crabbing  angle  has  been  estimated  (0.7  degrees  for  this  dataset),  the 
coordinates  of  the  phase  centres  can  be  derived:  for  each  ping,  the  array  is  tilted  w  ith 
the  crabbing  angle  as  shown  in  Figure  4. 13  and  the  derived  calculated  midpoints 
coordinates  will  then  be  used  for  the  motion  compensation. 
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x  1  q  3  Channels  and  Midpoints  coordinates 


Azimuth  [m] 


Figure  4. 13  Tilted  array:  in  blue,  the  calculated  element  positions  and  in  red  the  corresponding  midpoint 
positions  Note  that  scales  in  x-  and  y-direction  are  very  different, 


Using  the  derived  positions,  their  range  offsets  are  used  to  calculate  time  shifts  which 
are  then  applied  to  the  corresponding  element  signals.  The  result  of  this  time  shift  w  ill 
be  that  all  motion  corrected  midpoints  lie  along  the  same  line.  Once  this  motion 
correction  has  been  done,  we  can  display  the  common  midpoint  gathers  to  see  the 
improvement  as  compared  with  Figure  4.8: 


CMP  gather  0501 


CMP  gather  0601 


12  3  4  stack 


(b) 

Figure  4  14  Common  midpoints  gather,  (a)  is  without  motion  compensation,  (b)  is  with  motion 

compensation  We  can  notice  that  the  small  phase  shift  that  was  observed  without  motion 
compensation  is  now  much  smaller:  on  (b)  the  four  left-hand  curves  are  in  phase  and  the  stack 
curve  has  about  the  same  amplitude  as  the  four  others. 
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The  small  phase  shift  that  appeared  in  Figure  4.8  is  now  much  smaller.  As  for  the 
amplitude  of  the  fifth  curve,  it  now  equals  the  amplitude  of  the  first  four  curves. 

It  should  be  noted  here  that  the  rail  configuration  (i.e.  assuming  the  sway  was 
negligible)  allowed  us  to  derive  the  crabbing  angle.  On  a  more  conventional  platform 
(ship,  submarine,  AUV  ...)  the  DPCA  algorithm  needs  to  be  generalized  (including 
subsequent  midpoint/phase  centre  positions)  to  derive  the  crabbing  angle  and  sway 
simultaneously. 


4.3  Results  on  DERA/GESMA  rail  data  set  011004*  (75  m  range) 

These  data  image  an  area  with  both  a  sphere  and  a  cylinder  at  a  range  of  approximately 
75  m.  Figure  4. 1 5  shows  the  SAS  image  constructed  from  the  data  acquired  by  array  1 . 


70  72  74  76  78  80  82  84  86  88  90 

Range  [m] 

Figure  4  15  SAS  image  for  DERA/GLSMA  rail  data  01 1004,  array  1 


The  acquisition  parameters  are  listed  in  Table  4.4.  Information  on  the  origin  time  of  the 
data  and  the  height  of  the  platform  was  not  available.  In  order  to  determine  the  range 
bracket,  we  repeated  the  Stolt  migration  for  various  target  ranges,  optimising  the 
contrast  of  the  spherical  target.  An  estimate  of  the  platform  height  is  inferred  from  the 
shadow  length.  This  estimate  is  further  tuned  to  locate  the  seafloor  at  z  =  0  m  depth. 


Processing:  Software  id=  *DLRA_GESMA_raildata;  01 1004 
Interferometry:  Software  id  -  401 1004  20080 128’. 


3 
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Table  4  4  Acquisition  parameters  for  DERA/GESMA  rail  data  011004. 


Pulse  centre  frequency  fc 

150  kHz 

Pulse  bandwidth  b 

60  kHz 

Pulse  type 

Linear  frequency  modulated  pulse  (LFM) 

Pulse  duration 

4.0  ms 

Length  of  synthetic  aperture 

10.7  m 

Number  of  pings 

320 

Number  of  elements/ping 

64  (32  in  each  array) 

Separation  of  elements 

0  00834  m 

Displacement  between  pings 

0.0334  m 

Platform  height  relative  to  seafloor  at  target 

Estimated  7.5  m  (±  0  5  m) 

Platform  y-distance  to  targets 

Approximately  75  m  (±  0.2  m) 

Number  of  arrays 

2 

Baseline  vector  B  before  tilting 

(0,  0,  0.217)  m 

Baseline  tilt  angle 

6  degrees 

Origin  of  range  vector  r0 

Estimated  68  30  m  (±  0.2  m) 

In  order  to  assess  the  success  of  the  geometry-based  co-registration,  we  focus  on  a  part 
of  the  cylinder.  Figure  4.16  shows  that  the  geometrical  displacement  between  the 
master  and  co-registered  slave  image  is  small,  and  thus  indicates  the  success  of  the 
geometry-based  co-registration. 


HnW  taugn  (dB)  iton  In#  |4< 


Figure  4  16  Detail  of  the  SAS  image  of  the  cylinder  illustrated  in  Figure  4  1 5,  as  used  for  interferometry. 

Shown  are  the  master  image  obtained  with  array  1  on  the  left  and  the  co-registered  (slave) 
image  from  array  2  on  the  nght 


As  suggested  by  Figure  4. 16,  the  co-registered  master  and  slave  images  are  highly 
correlated.  This  is  confirmed  by  the  correlation  map  (Figure  4.17).  A  correlation  value 
equal  to  one  indicates  that  the  master  and  slave  images  are  perfectly  correlated 
Low  correlation  values  are  observed,  as  expected,  in  the  shadow  regions  (e.g.  behind 
the  sphere).  Multipath  propagation  is  another  factor  causing  data  to  decorrelate. 

The  correlation  map  is  an  important  measure  for  quality  control.  If  data  are  poorly 
correlated,  phase  estimates  are  likely  to  be  uncertain.  Then,  the  height  estimates  become 
unreliable. 
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Figure  4,17  Correlation  map,  indicating  local  correlation  coefficients  between  the  master  image  and  the  co¬ 
registered  slave  image.  High  correlation  is  required  to  achieve  good  height  estimates. 


Figure  4. 18  shows  the  interferometric  height  estimates.  Both  the  sphere  and  the  cylinder 
can  be  easily  recognised.  Height  values  with  a  reflectivity'  of  -45  dB  relative  to  the 
maximum  reflectivity  have  been  set  to  zero. 


Height  [m] 


Figure  4.18  Height  estimate  derived  from  the  interferogram,  showing  both  the  sphere  and  the  cylinder. 


Figure  4.19  and  Figure  4.20  focus  on  the  sphere  and  the  cylinder  and  show  2D  and  3D 
views  of  reflectivity.  The  image  of  the  sphere  has  two  highlights.  The  first  highlight 
corresponds  to  the  direct  path,  whereas  the  second  highlight  is  attributed  to  the  seabed 
multipath.  The  travel  time  of  this  multipath  corresponds  to  an  imaginary  target  located 
below  the  seabed.  This  imaginary  target  is  a  mirror  image  of  the  sphere.  The  3D 
reflectivity  map  shows  that  interferometry  indeed  locates  the  multipath  highlight  below 
the  seabed,  and  is  therefore  indicative  for  the  resolution  of  the  height  estimates. 
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Figure  4. 1 9  2D  (top)  and  3D  (bottom)  reflectivity  images  of  the  sphere 


In  general,  multipath  propagation  is  a  problem  for  interferometry.  Reliable  height 
estimates  cannot  be  obtained  when  a  pixel  has  a  direct  and  a  multipath  scattering 
contribution  of  comparable  strength.  We  expect  that  the  noise  in  the  height  estimates  for 
both  the  sphere  and  the  cylinder  is  at  least  partly  caused  by  multipath. 

Multipath  propagation  is  especially  relevant  at  low  grazing  angles. 
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Figure  4  20  2D  (top)  and  3D  (bottom)  reflectivity  images  of  the  cylinder 
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4.4  DERA/GESMA  rail  data  set  051003*  (25  m  range) 

This  data  set  images  an  area  with  a  target  at  a  range  of  approximately  25  m.  The  target 
is  a  sphere  w  ith  a  diameter  of  approximately  1.0  m.  As  a  result  of  the  short  range  in 
comparison  with  the  01 1004  data  set,  the  grazing  angle  is  larger.  We  expect  therefore 
that  multipath  will  be  less  significant.  Figure  4.21  shows  the  SAS  image  from  array  1 . 


Range  [m] 

Figure  4.21  SAS  image  from  DERA/GESMA  rail  data  051003. 

The  acquisition  parameters  are  listed  in  Table  4.5.  Again,  the  range  is  estimated  by 
optimizing  the  contrast  in  the  image.  The  height  of  the  platform  is  obtained  using  the 
shadow  and  fine-tuned  to  locate  the  background  at  z  =  0. 


Table  4.5  Acquisition  parameters  for  DERA/GESMA  rail  data  051003. 


Pulse  centre  frequency  fc 

150  kHz 

Pulse  bandwidth  b 

60  kHz 

Pulse  type 

Linear  frequency  modulated  pulse  (LFM) 

Pulse  duration 

4.0  ms 

Length  of  synthetic  aperture 

10.7  m 

Number  of  pings 

320 

Number  of  elements/ping 

64  (32  in  each  array) 

Separation  of  elements 

0.00834  m 

Displacement  between  pings 

0.0334  m 

Platform  height  relative  to  seafloor  at  target 

Estimated  9.2  m  (±  0.5  m) 

Platform  y-distance  to  targets 

Approximately  25  m  (±  0  2  m) 

4  Processing.  Software  id  =  ‘DERAGESMAraildata,  05 1 003. 
Interferometry:  Software  id  =  *051003  20080128’. 
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Table  4.5  Acquisition  parameters  for  DERA/GESMA  rail  data  051003  (continued). 


Number  of  arrays 

2 

Baseline  vector  B  before  tilting 

(0,  0,  0.217)  m 

Baseline  tilt  angle 

14  degrees 

Origin  of  range  vector  ro 

Estimated  13  50  m  (±  0.2  m) 

Details  of  the  background  imaged  in  both  the  master  and  slave  image  show  that  these 
are  highly  correlated  (Figure  4.22),  and  this  is  confirmed  by  the  correlation  map 
(Figure  4.23).  The  data  are  highly  correlated  except  on  the  target  itself. 


Figure  4.22  Detail  of  the  background  in  the  SAS  image  illustrated  in  Figure  4  2 1  Shown  are  the  master 
image  from  array  1  and  the  co-registered  image  from  array  2. 


Correlation  map 


Azimuth  [m] 


Figure  4.23  Correlation  map,  indicating  local  correlation  coefficients  between  the  master  image  and  the 

co-registered  slave  image  of  the  sphere  and  a  patch  of  surrounding  seafloor  High  correlation  is 
required  to  achieve  good  height  estimates. 


The  interferometric  height  estimates  (Figure  4.24  and  Figure  4.25)  show  that  the 
background  is  very  flat  and  well  resolved.  On  the  other  hand,  the  height  of  the  target  is 
not  well  resolved.  This  is  presumably  related  to  phase  unwrapping  problems. 

These  unwrapping  errors  may  be  caused  by  the  discontinuity  in  the  height  at  the  edge  of 
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the  sphere.  Assuming  that  the  subsidence  of  the  target  can  be  neglected,  the  height  of 
the  normal  incidence  target  highlight  is  located  approximately  0.6  m  above  the  seafloor, 
about  10  cm  higher  than  the  sphere  midpoint  due  to  the  incidence  angle. 

The  corresponding  path  differences  to  array  1  and  2  are  of  similar  order  of  magnitude  as 
the  wavelength  (0.01  m).  Therefore,  we  expect  that  the  negative  height  values  are  phase 
unwrapping  artefacts.  Another  possibility  is  multipath  propagation  as  discussed  in  the 
previous  section. 

Height  [m] 


Azimuth  [m] 


Figure  4.24  Height  estimate  derived  from  the  interferogram  (after  phase  unwrapping). 
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Figure  4.25  2D  (a)  and  3D  (b)  reflectivity  images  of  the  sphere 
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5  Conclusions  and  outlook 


The  preceeding  chapters  have  described  the  status  of  interferometric  SAS  processing  at 
TNO,  The  current  implementation  uses  SAS  processing  as  an  initial  step.  Then, 
interferometric  processing  is  applied  to  the  SAS  images  to  infer  the  height.  Experiments 
with  both  simulated  and  measured  data  provided  promising  and  interesting  results. 

These  gave  insight  into  the  prospects  -  such  as  valid  height  values  corresponding  to  the 
targets  and  for  the  bathymetry  -  but  also  helped  us  to  identify  limitations  of  SAS 
interferometry. 

The  most  important  limitation  is  the  variability  of  the  height  estimates  in  both  the 
simulated  and  measured  data.  This  addresses  a  fundamental  issue,  i.e.  that  the  SAS 
processing  needs  to  be  phase-preserving .  In  the  presence  of  strong  scatterers,  imaging 
artefacts  ( sidelobes )  in  the  vicinity  of  the  scatterer  cause  a  considerable  amount  of  noise 
in  the  phase  information.  Furthermore,  interferometry  is  very  sensitive  to  platform 
motion  errors,  and  height  estimates  are  also  unreliable  in  the  case  of  multipath 
propagation  (height  ambiguity  with  range). 

Further  development  of  SAS  interferometry'  should  be  aimed  at  improving  the  reliability' 
of  the  phase  information.  The  following  aspects  are  aimed  at  resolving  the  above- 
mentioned  bottlenecks: 

•  We  need  to  consider  the  interferometric  processing.  Instead  of  starting  with  the  SAS 
imaging,  which  aims  at  optimising  the  contrast  in  the  images,  specific 
interferometric  imaging  needs  to  be  considered.  This  processing  should  include  a 
compensation  for  highlights  to  suppress  sidelobes. 

•  In  addition,  the  higher  sensitivity'  to  motion  errors  requires  to  reconsider  the 
optimum  integration  angle  (or  synthetic  aperture  length).  A  shorter  aperture 
decreases  the  sensitivity  to  motion  errors.  This  should  be  balanced  with  the  loss  in 
azimuth  resolution.  The  results  obtained  in  different  sub-apertures  can  be  combined, 
to  improve  the  S  ^N  ratio  of  the  phase  estimates. 

•  It  has  been  shown  that  the  accuracy  of  the  instantaneous  phase  estimate  decreases 
with  increasing  bandwidth  (Figure  2.4).  As  a  consequence,  there  is  a  tradeoff 
between  the  accuracy  of  the  interferometric  height  estimate  and  the  range  resolution 
that  can  be  achieved,  which  advocates  processing  at  reduced  bandwidth.  In  general, 
the  interferometric  resolution  will  be  significantly  lower  than  in  the  reflectivity 
images. 

•  Finally,  a  more  sophisticated  method  for  phase  unwrapping  has  to  be  implemented. 
Phase  unwrapping  errors  lead  to  erroneous  height  estimates.  Special  attention  is 
required  to  resolve  these  ambiguities  when  there  are  discontinuities  in  height 

(e.g.  when  a  spherical  object  is  located  on  the  seabed). 

A  remaining  challenge  is  to  obtain  reliable  phase  information  in  the  presence  of 
multipath  propagation,  which  is  also  referred  to  as  overlay. 

The  developed  algorithms  will  eventually  be  applied  to  low-frequency  buried-mine 
data.  It  is  expected  that  distinguishing  bottom  echoes  from  echoes  from  below  will 
improve  detection  performance.  An  experiment  on  objects  buried  in  mud  will  most 
likely  take  place  in  2009.  Adjustments  and  modifications  to  the  software  will  no  doubt 
be  necessary,  but  the  interferometric  toolbox  is  now  available  and  operational. 


TNO  report  |  TNO-DV  2008  A176 


42/44 


6  References 


[1]  Rosen,  P.A.,  Hensley,  S.  Joughin,  I.R.,  Li,  F.K.,  Madsen,  S.N.,  Rodriguez,  E.,  and 
Goldstein,  R.M.,  “ Synthetic  aperture  radar  interferometry". 

Proceedings  of  the  IEEE,  88(3),  333-382,  2000. 

[2]  Bamler,  R.,  and  Hartl,  P.,  “ Synthetic  aperture  radar  interferometry" , 

Inverse  Problems  14,  RI-R54,  1998. 

[3]  Griffiths,  H.D.,  et  al., 44 Interferometric  synthetic  aperture  sonar  for  high- 
resolution  3-D  mapping  of  the  seabed ", 

IEE  Proc.  Radar  Sonar  Navigation  144(2):  96-103,  1997. 

[4]  Saebo,  T.O.,  Hansen,  R.E.,  Callow,  H.J.,  “ Height  estimation  on  wideband 
synthetic  aperture  sonar:  experimental  results  from  InSAS-2000" ,  Oceans  2005. 

[5]  Bonifant,  W.W.,  “ Interferometric  synthetic  aperture  sonar  processing" , 

Msc.  Thesis,  Georgia  Institute  of  Technology,  1999. 

[6]  Barclay,  J.P.,  44 Interferometric  synthetic  aperture  sonar  design  and  performance" , 
PhD  Thesis,  University  of  Canterbury  ,  2006. 

[7]  Hansen,  R.E.,  Saebo,  T.O.,  Callow,  II. J.  and  Hagen,  P.E.,  “ Synthetic  aperture 
sonar  processing  for  the  HUGIN  A  UV\  Oceans  2005. 

[8]  Hagen,  P.E.,  Hansen,  R.E.,  and  Langli,  B.,  44 Interferometric  synthetic  aperture 
sonar  for  the  HUGIN  1000-MRAUV',  UDT  Pacific  2006. 

[9]  Vernier,  A.,  “ Interferometric  synthetic  aperture  sonar  for  the  detection  of  buried 
objects"  Msc.  Thesis,  ENSIETA,  2007. 

[10]  Barnes,  A.  E.,  “A  tutorial  on  complex  trace  analysis ", 

Geophysics,  72,  pp.  W33-W43,  2007. 

[11]  Sheriff,  R.E.,  44 Encyclopedic  dictionary'  of  applied  geophysics,  4th  ed\ 

Society  for  Exploration  Geophysicists,  2002. 

[12]  Goldstein,  R.M.,  Zebker,  H.A.,  and  Werner,  C.L., 44 Satellite  radar  interferometry: 
two-dimensional  phase  unwrapping" , 

Radio  Science,  23,  713-720,  1988. 

[13]  Ghighlia,  D.C.,  and  Pritt,  M.D., 44 Two-dimensional  phase  unwrapping:  Theory , 
Algorithms,  and  Software ",  Wiley,  New  York,  1998. 

[14]  Chen,  C.W.,  and  Zebker,  H.A., 44 Network  approaches  to  two-dimensional  phase 
unwrapping" ,  Journal  of  the  Optical  Society  of  America  A,  17(3),  401-414,  2000. 

[15]  Chen,  C.W.,  and  Zebker,  H.A., 44 Two-dimensional  phase  unwrapping  with  use  of 
statistical  models  for  cost  functions  in  nonlinear  optimization ", 

Journal  of  the  Optical  Society  of  America  A,  18(2),  338-351,  2001. 

[  1 6]  Chen,  C.W.,  and  Zebker,  H.A., 44 Phase  unwrapping  for  large  SAR 

Interferograms:  Statistical  segmentation  and  generalized  netw  ork  models ", 

IEEE  Transactions  on  Geoscience  and  Remote  Sensing,  40,  1709-1719. 

[17]  Groen,  J., 44 Adaptive  motion  compensation  in  sonar  array  processing", 

PhD.  Thesis  TU  Delft,  2006. 

[18]  Stolt,  R.H., 44 Migration  by  Fourier  transform".  Geophysics,  43,  23-48,  1978. 


TNO  report  |  TNO-DV  2008  A176 


43/44 


[19]  Yilmaz,  O.,  “ Seismic  data  analysis,  vol  r\ 

Society  for  Exploration  Geophysicists,  2002. 

[20]  Belletini,  A.  and  Pinto,  M.  A.,  “ Theoretical  accuracy  of  synthetic  aperture  sonar 
micronavigation  using  a  displaced  phase-center  antenna  ", 

IEEE  Journal  of  Oceanic  Engineering,  27,  780-789,  2002. 


TNO  report  |  TNO-DV  2008  A176 


44/44 


7  Signature 


The  Hague,  October  2008  TNO  Defence,  Security  and  Safety 


ONGERUBRICEERD 


REPORT  DOCUMENTATION  PAGE 

(MOD-NL) 


1  DEFENCE  REPORT  NO  (MOD-NL) 

TD2008-0067 

2  RECIPIENT'S  ACCESSION  NO 

3  PERFORMING  ORGANIZATION  REPORT  NO 

TNO-DV  2008  A 176 

4  PROJECT/TASKWORK  UNIT  NO 

015.34736 

5.  CONTRACT  NO 

6.  REPORT  DATE 

October  2008 

7.  NUMBER  OF  PAGES 

44  (excl  RDP  &  distribution  list) 

8  NUMBER  OF  REFERENCES 

20 

9  TYPE  OF  REPORT  AND  DATES  COVERED 

Final 

10.  TITLE  AND  SUBTITLE 

The  prospects  of  SAS  interferometry  for  detection  and  classification 

11.  AUTHOR(S) 

Dr  R.  van  Vossen 
B.A.J.  Quesson 
Dr  J.C.  Sabel 

12.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

TNO  Defence,  Security  and  Safety,  P.O.  Box  96864,  2509  JG  The  Hague,  The  Netherlands 
Oude  Waalsdorperweg  63,  The  Hague,  The  Netherlands 

13  SPONSORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Royal  Netherlands  Army 

14  SUPPLEMENTARY  NOTES 

The  classification  designation  Ongerubriceerd  is  equivalent  to  Unclassified,  Stg.  Confidentieel  is  equivalent  to 
Confidential  and  Stg.  Geheim  is  equivalent  to  Secret. 

15.  ABSTRACT  (MAXIMUM  200  WORDS  (1044  BYTE)) 

Processing  for  interferometric  synthetic  aperture  sonar  has  been  developed,  and  tested  with  both  measured  and 
simulated  data.  Interferometry  is  based  on  data  from  two  vertically  separated  receive  arrays.  Subtle  phase  differences 
between  the  images  from  both  arrays  provide  information  on  the  relative  height  of  objects  in  the  observed  scene.  The 
results  give  confidence  in  the  processing  and  provide  insight  into  the  limitations  and  options  for  improvement  of  the 
current  processing  suite.  Eventually,  the  developed  processing  is  expected  to  improve  the  detection  of  objects  buried 


in  the  sea  bottom. 

16  DESCRIPTORS 

IDENTIFIERS 

Synthetic  aperture  sonar,  interferometry,  buried  mines 

17a  SECURITY  CLASSIFICATION  17b  SECURITY  CLASSIFICATION 

(OF  REPORT)  (OF  PAGE) 

Ongerubriceerd  Ongerubriceerd 

17c.  SECURITY  CLASSIFICATION 

(OF  ABSTRACT) 

Ongerubriceerd 

18.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Unlimited  Distribution 

17d  SECURITY  CLASSIFICATION 

(OF  TITLES) 

Ongerubriceerd 

ONGERUBRICEERD 


Distribution  list 


The  following  agencies/people  will  receive  a  complete  copy  of  the  report. 

I  DMO/SC-DR&D 

standaard  inclusief  digitale  versie  bijgeleverd  op  cd-rom 

2/3  DMO/DR&DICennistransfer 

4  Programma-'projectbegeleider  Defensie 
LTZ1  drs.  R.P.A.  Dekeling 

5  DMO/DWS&B  RZS&B  Sensor-  en  Wapentechnologie  (SWT) 
LTZl  B.E.A.  Kerstens,  MSc 

6/8  Bibliotheek  KMA 

9  Programmaleider  TNO  Defensie  en  Veiligheid,  dr.ir.  J.C.  Sabel 

10/11  TNO  Defensie  en  Veiligheid,  vestiging  Den  Haag, 

Archief 

12/17  TNO  Defensie  en  Veiligheid,  vestiging  Den  Haag, 

Business  Unit  Waamemingssystemen, 

dr.  R.  van  Vossen 

B.A.J.  Quesson 

ir.  F.P.G.  Driessen 

dr.ir.  G.  Blacquiere 

dr.  S.P.  Beerens 

dr.ir.  A.L.D.  Beckers 


The  following  agencies/people  will  receive  the  management  summary  and  the 
distribution  list  of  the  report. 


4  ex. 

DMO'SC-DR&D 

1  ex. 

DMO^ressort  Zeesystemen 

1  ex. 

DMO/ressort  Landsystemen 

1  ex. 

DMO/ressort  Luchtsystemen 

2  ex. 

BS/DS  DOBBP/SCOB 

1  ex. 

MIVD  AAR/BMT 

1  ex. 

StafCZSK 

1  ex. 

StafCLAS 

1  ex. 

StafCLSK 

1  ex. 

StafKMar 

1  ex. 

TNO  Defensie  en  Veiligheid,  Algemeen  Directeur, 
ing.  J.V.  Elsendoom 

1  ex. 

TNO  Defensie  en  Veiligheid,  Directie 

Directeur  Operaties,  ir.  C.  Eberwijn 

1  ex. 

TNO  Defensie  en  Veiligheid.  Directie 

Directeur  Kennis,  prof  dr.  P.  Werkhoven 

1  ex. 

TNO  Defensie  en  Veiligheid,  Directie 

Directeur  Markt,  G.D.  Klein  Baltink 

1  ex. 

TNO  Defensie  en  Veiligheid,  vestiging  Den  Haag, 

Manager  Waamemingssystemen  (operaties),  ir.  B.  Dunnebier  PDeng 

1  ex. 

TNO  Defensie  en  Veiligheid,  vestiging  Den  Haag, 

Manager  Informatie  en  Operaties  (operaties),  ir.  P.  Schulein 

1  ex. 

TNO  Defensie  en  Veiligheid,  vestiging  Rijswijk,  daama  reserve 

Manager  Bescherming,  Munitie  en  Wapens  (operaties),  ir.  P.J.M.  Elands 

1  ex. 

TNO  Defensie  en  Veiligheid,  vestiging  Rijswijk, 

Manager  BC  Bescherming  (operaties),  ir.  R.J. A.  Kersten 

1  ex. 

TNO  Defensie  en  Veiligheid,  vestiging  Soesterberg, 

Manager  Human  Factors  (operaties),  drs.  H.J.  Vink 

